function [T, cu, cl, r, V_u, rho_u, P_l, cr_topo] = dynamicPlume(mu_l, mu_u, rho_l, rho_lm, delTemp, V_flux, t_f)
%% Mantle Plume Flow Model, 2023 (Molitor, Royden, and Jagoutz)
%Created: January, 2023
%Author: Zachary Molitor
%Dependencies: bandthree.m

%January, 4, 2023, fixed volume issue originating on line 101, V_flux/(r_0)
%to V_flux/(2*pi*r_0)

%January, 20, 2023, Wiki (Leigh) fixed a term in equation 4 (rho_u replaced
%with rho_l) and also incoporated water loading into the main isostatic
%balance of the code

%Description: Taking input viscosities, densities, volume flux, run time,
%and a lower mantle flux term, this function calculates the dynamic
%topography of a plume. The plume conduit is already in place and has
%reached the base of the lithosphere by t = 0 of the model. After t = 0,
%plume material is fluxed from the radially symmetric, cylindrical plume conduit into
%the surrounding upper mantle, creating a plume channel between the
%lithosphere and the rest of the upper mantle. The lithosphere and the interface
%between the upper and lower mantle are deflected to accomodate the new
%channel layer being injected, and also to accomodate the pressure gradient
%that drives the lateral flow away from the plume conduit.

%INPUTS:
%mu_l: Viscosity of the non-plume upper mantle (Pa*s)
%mu_u: Viscosity of the plume material (Pa*s)

%rho_l: Density of average upper mantle (kg/m^3)
%rho_lm: Density of the lower mantle, below about 200-300km depth (kg/m^3)
%delTemp: The difference in temperature between the plume material and the
%   ambient upper mantle temperature on a global scale. 
%   This is used to calculate the density of the plume material.(C/K)

%V_flux: Volume flux of the plume. (km^3/Myr)

%t_f: Run time of the model, generally in millions of years, time step is
%   10^10s, so needs to be signficantly larger than the time step. (yrs)

%f_ltolm: Fractional flux from the upper mantle into the lower mantle
%   representing loss of material as the mantle root grows. A value from 0 to
%   1, where 0 is no loss and 1 is total loss of any mantle root material.


%OUTPUTS:
%T: The calculated topography from the model, based on inputs. (m)
%cu: The thickness of the plume channel in the upper mantle (m)
%cl: The thickness of the lower upper mantle layer, where cl+cu is the
%total thickness of the upper mantle (m)
%r: distance from mantle plume center array (m)


%% Set constants and initial variables
%Variables of note: mu's, rho's, delTemp, volume flux

addpath('C:\Users\moza0\Dropbox (MIT)\plumeflowmodel\addaxis6')
% gravitational const (m/s2)
g = 9.81;

% density of water
rhow = 1000;

% coefficient of thermal expansion (1/Kelvin)
% alpha = 3e-5; %alternate value
alpha = 2e-5; %changed 1/30/2023, likely more rep of lithospheric mantle

% radial grid spacing
delr = 5; %kilometers
% delr = 10; %kilometers

%initial channel thickness (realistically a very small number or 0) (m)
cu_0 = 0;

% Plume head total radius (kilometers)
r_o = 100; %comp to hongze
% r_o = 100; %rs, iceland
% r_o = 400; %testing, added 1/30/23

%total distance of radial array
r_f = 5000;

% time step
delt = 1e10; %seconds
% delt = 5e8; %seconds

%seconds per year
speryr = 3.154e7;

%negative fractional flux term, represents loss of upper mantle to lower
%mantle as a function of the size of the upper mantle root
% f_ltolm = 0; %Choose value from 0 to 1, 0 is no loss, 1 is total loss

%% Convert inputs to SI units and calculate any additional starting variables
% IMPORTANT: time variable below, set total run time of model

% initial distance to base of upper mantle (convert to meters)
% cl_0 = cl_0*1000;
% Plume head total radius (convert to meters)
r_o = r_o*1000;
r_f = r_f*1000;
% Total volume flux V_o*ro (km^3/Myr to m^3/s)
V_flux = V_flux * 3.154e-5;
% radial grid spacing to meters
delr = delr * 1000;

%density contrast
delRho = delTemp*rho_l*alpha;

% Flux at plume radius
V_o = V_flux./(2*pi*r_o); %deprecated with Q

% density within channel
rho_u = rho_l - delRho;

% radial array (in meters)
r = delr:delr:r_f; %changed for Q, now starts model at 0 km instead of plume radius


% initial thickness of horizontal plume channel (meters)
cu = cu_0*ones(size(r));
% setting up array for temp storage of new channel thickness, used in each
% time step
cu_new = zeros(size(r));
cl_new = zeros(size(r));

%set initial thickness of asthenosphere layer (one layer in simple case, but changes in thickness still happen)
% cl = cl_0*ones(size(r));

%initial topography over grid
T = zeros(size(r));

t_0 = 0; % start time
% t_f = speryr*2.7e7; % final time to finish loop,iceland
% t_f = speryr*1.5e7; % for time series (iceland params)
% t_f = speryr*2.9e7; %rs
% t_f = speryr*2e7; %azores
% t_f = speryr*2e7; %testing
t_f = speryr*t_f;

% arbitrary constants for calculation
w = delt./((delr.^2).*2);

a_0 = 30;
a_0 = a_0*1000;
a_0_array = a_0*ones(size(r));
cl = a_0_array;

V_u = zeros(size((t_0):delt:(t_f)));
Q = V_flux./(pi.*(r_o.^2));
% Q = V_flux./(pi.*(r.^2));
italic_f = 1e-18;
P_l = 0;

%% Read and interpolate crustal thickness data

D = readtable("IcelandCrustalThickness\Iceland_crustalthickness_weiretal2001.xlsx");
D = table2array(D);
r_cr = D(:,1).*1000;
cr_topo = D(:,5);
cr_topo = interp1(r_cr,cr_topo,r,'pchip');
% cr_topo = smooth(cr_topo,300);
% cr_topo = transpose(cr_topo);

%% Implicit finite difference
% bandthree.m used to find matrix solution

for t = 0:delt:(t_f)

    Told = T; %temp save for old topography for dT calculation
    
    % Arbitrary constants for upper velocity flux and calculations
    % below. (viscosity dependent on channel thickness also maybe have check for channel thickness > 50 and make it super strong)
%     mu_u_scaling = mu_u.*(10./(cu+.0001)); %scale plume viscosity with channel thickness, thicker channels have relatively lower viscosity due to slower cooling

    % CONSTANT MU_U
    B_u = (((rho_u-rhow).*g.*((mu_l*(cu.^4)) + (4.*mu_u.*cl.*(cu.^3)))) + ((rho_l-rhow).*g.*3.*mu_u.*(cl.^2).*(cu.^2)))...
        ./(12.*mu_u.*(cu.*mu_l+mu_u.*cl)); 
    D_u = (rho_u-rho_l).*g.*((cl.^2).*(cu.^2))...
        ./(4.*(cu.*mu_l+mu_u.*cl)); 

    B_l = (((rhow-rho_l).*g.*((-1.*mu_u.*(cl.^4)) + (-4.*mu_l.*cu.*(cl.^3)))) - ((rhow-rho_u).*g.*3.*mu_l.*(cl.^2).*(cu.^2)))...
        ./(12.*mu_l.*(cu.*mu_l+mu_u.*cl)); 
    D_l = (rho_l-rho_u).*g.*(((-1.*mu_u.*(cl.^4)+(-4.*mu_l.*cu.*(cl.^3))))...
        ./(12.*mu_l.*(cu.*mu_l+mu_u.*cl))); 

%     % SCALING MU_U
%     B_u = (((rho_u-rhow).*g.*((mu_l*(cu.^4)) + (4.*mu_u_scaling.*cl.*(cu.^3)))) + ((rho_l-rhow).*g.*3.*mu_u_scaling.*(cl.^2).*(cu.^2)))...
%         ./(12.*mu_u_scaling.*(cu.*mu_l+mu_u_scaling.*cl)); 
%     D_u = (rho_u-rho_l).*g.*((cl.^2).*(cu.^2))...
%         ./(4.*(cu.*mu_l+mu_u_scaling.*cl)); 
%     
%     B_l = (((rhow-rho_l).*g.*((-1.*mu_u_scaling.*(cl.^4)) + (-4.*mu_l.*cu.*(cl.^3)))) - ((rhow-rho_u).*g.*3.*mu_l.*(cl.^2).*(cu.^2)))...
%         ./(12.*mu_l.*(cu.*mu_l+mu_u_scaling.*cl)); 
%     D_l = (rho_l-rho_u).*g.*(((-1.*mu_u_scaling.*(cl.^4)+(-4.*mu_l.*cu.*(cl.^3))))...
%         ./(12.*mu_l.*(cu.*mu_l+mu_u_scaling.*cl))); 

    
    
    %Arbitrary variables for topography finite difference
    B_t = B_u+B_l;
    D_t = D_u+D_l;

%     P_l = ((rho_l-rhow).*g.*T) + (rho_u-rho_l).*g.*cu;
    %INPUT FLUX BC
    %initialize variables for band matrix
    A = zeros(size(B_u));
    B = zeros(size(B_u));
    C = zeros(size(B_u));
    f = zeros(size(B_u));


    %Q term
    %(Q*(r_o.^2./100000.^2).*exp(-r(1).^2./100000.^2))
    
    %Conditional check for cl thickness (for model including E)
    if (cl(1) < 1) %dT/dt = dcu/dt-EPl
        f(1) = T(1) + (w).*(1/r(1)).*(r(2).*D_u(2)+r(1).*D_u(1)).*(cu(2)-cu(1))...
            + (delt.*(Q*(r_o.^2./20000.^2).*exp(-r(1).^2./20000.^2))-(delt.*italic_f.*(rho_u-rho_l).*g.*cu(1)));
        
        B(1) = 1+(delt.*italic_f.*(rho_l-rhow).*g)+(w).*(1/r(1)).*(r(2).*B_u(2)+r(1).*B_u(1));
        C(1) = (-w).*(1/r(1)).*(r(2).*B_u(2)+r(1).*B_u(1));
    else %standard dT/dt = dcu/dt+dcl/dt
        f(1) = T(1) + (w).*(1/r(1)).*(r(2).*D_t(2)+r(1).*D_t(1)).*(cu(2)-cu(1))...
            + (delt.*(Q*(r_o.^2./20000.^2).*exp(-r(1).^2./20000.^2))-(delt.*italic_f.*(rho_u-rho_l).*g.*cu(1)));

        B(1) = 1+(delt.*italic_f.*(rho_l-rhow).*g)+(w).*(1/r(1)).*(r(2).*B_t(2)+r(1).*B_t(1));
        C(1) = (-w).*(1/r(1)).*(r(2).*B_t(2)+r(1).*B_t(1));
    end

    %Assign endmember values of band matrix (unaffected by cl value)
    A(1) = 0;
    A(size(r,2)) = 0;
    C(size(r,2)) = 0;
    B(size(r,2)) = 1;
    f(size(r,2)) = 0;



    % create center and upper bands
    for i = 2:(size(r,2)-1)
        if (cl(i) < 1) %dT/dt = dcu/dt-EPl
            C(i) = (-w).*(1/r(i)).*(r(i+1).*B_u(i+1) + r(i).*B_u(i));
            B(i) = 1 + (delt.*italic_f.*(rho_l-rhow).*g)+(w).*(1/r(i)).*(r(i+1).*B_u(i+1)+2.*r(i).*B_u(i)+r(i-1).*B_u(i-1));
            A(i) = (-w).*(1/r(i)).*(r(i).*B_u(i)+r(i-1).*B_u(i-1));
            if (r(i) < r_o)
                f(i) = T(i) + ...
                    cu(i+1).*(w).*(1/r(i)).*(r(i+1).*D_u(i+1) + r(i).*D_u(i))...
                    + cu(i).*(-w).*(1/r(i)).*(r(i+1).*D_u(i+1)+2.*r(i).*D_u(i)+r(i-1).*D_u(i-1))...
                    + cu(i-1).*(w).*(1/r(i)).*(r(i).*D_u(i)+r(i-1).*D_u(i-1))...
                    + (delt.*(Q*(r_o.^2./20000.^2).*exp(-r(i).^2./20000.^2))-(delt.*italic_f.*(rho_u-rho_l).*g.*cu(i)));
            else
                f(i) = T(i) + ...
                    cu(i+1).*(w).*(1/r(i)).*(r(i+1).*D_u(i+1) + r(i).*D_u(i))...
                    + cu(i).*(-w).*(1/r(i)).*(r(i+1).*D_u(i+1)+2.*r(i).*D_u(i)+r(i-1).*D_u(i-1))...
                    + cu(i-1).*(w).*(1/r(i)).*(r(i).*D_u(i)+r(i-1).*D_u(i-1))...
                    + delt.*(-(italic_f.*((rho_u-rho_l).*g.*cu(i))));
            end
        else %standard dT/dt = dcu/dt+dcl/dt
            C(i) = (-w).*(1/r(i)).*(r(i+1).*B_t(i+1) + r(i).*B_t(i));
            B(i) = 1 + (delt.*italic_f.*(rho_l-rhow).*g)+(w).*(1/r(i)).*(r(i+1).*B_t(i+1)+2.*r(i).*B_t(i)+r(i-1).*B_t(i-1));
            A(i) = (-w).*(1/r(i)).*(r(i).*B_t(i)+r(i-1).*B_t(i-1));
            if (r(i)<r_o)
                f(i) = T(i) + ...
                    cu(i+1).*(w).*(1/r(i)).*(r(i+1).*D_t(i+1) + r(i).*D_t(i))...
                    + cu(i).*(-w).*(1/r(i)).*(r(i+1).*D_t(i+1)+2.*r(i).*D_t(i)+r(i-1).*D_t(i-1))...
                    + cu(i-1).*(w).*(1/r(i)).*(r(i).*D_t(i)+r(i-1).*D_t(i-1))...
                    + (delt.*(Q*(r_o.^2./20000.^2).*exp(-r(i).^2./20000.^2))-(delt.*italic_f.*(rho_u-rho_l).*g.*cu(i)));
            else
                f(i) = T(i) + ...
                    cu(i+1).*(w).*(1/r(i)).*(r(i+1).*D_t(i+1) + r(i).*D_t(i))...
                    + cu(i).*(-w).*(1/r(i)).*(r(i+1).*D_t(i+1)+2.*r(i).*D_t(i)+r(i-1).*D_t(i-1))...
                    + cu(i-1).*(w).*(1/r(i)).*(r(i).*D_t(i)+r(i-1).*D_t(i-1))...
                    + delt.*(-(italic_f.*((rho_u-rho_l).*g.*cu(i))));
            end
        end
    end
    
    %calculate new topography using implicit finite matrix
    Tnew = bandthree(size(r,2),A,B,C,f);

    T = Tnew; %assign new values of T and use in calculating cu and cl

    % %IMPLICIT CU CALCULATION (try an explicit calculation)
    % Ac = zeros(size(B_u));
    % Bc = zeros(size(B_u));
    % Cc = zeros(size(B_u));
    % fc = zeros(size(B_u));
    % fc(1) = cu(1) + (w).*(1/r(1)).*(r(2).*B_u(2)+r(1).*B_u(1)).*(T(2)-T(1))...
    %     + (delt.*Q(1));
    % 
    % % setting boundaries of matrix below
    % Bc(1) = 1+(w).*(1/r(1)).*(r(2).*D_u(2)+r(1).*D_u(1));
    % Cc(1) = (-w).*(1/r(1)).*(r(2).*D_u(2)+r(1).*D_u(1));
    % Ac(1) = 0;
    % Ac(size(r,2)) = 0;
    % Cc(size(r,2)) = 0;
    % Bc(size(r,2)) = 1;
    % fc(size(r,2)) = 0;
    % 
    % % create center and upper bands
    % for i = 2:(size(r,2)-1)
    % 
    %     Cc(i) = (-w).*(1/r(i)).*(r(i+1).*D_u(i+1) + r(i).*D_u(i));
    %     Bc(i) = 1 + (w).*(1/r(i)).*(r(i+1).*D_u(i+1)+2.*r(i).*D_u(i)+r(i-1).*D_u(i-1));
    %     Ac(i) = (-w).*(1/r(i)).*(r(i).*D_u(i)+r(i-1).*D_u(i-1));
    %     if (r(i)<r_o)
    %         fc(i) = cu(i) + ...
    %             T(i+1).*(w).*(1/r(i)).*(r(i+1).*B_u(i+1) + r(i).*B_u(i))...
    %             + T(i).*(-w).*(1/r(i)).*(r(i+1).*B_u(i+1)+2.*r(i).*B_u(i)+r(i-1).*B_u(i-1))...
    %             + T(i-1).*(w).*(1/r(i)).*(r(i).*B_u(i)+r(i-1).*B_u(i-1))...
    %             + (delt.*Q);
    %     else
    %         fc(i) = cu(i) + ...
    %             T(i+1).*(w).*(1/r(i)).*(r(i+1).*B_u(i+1) + r(i).*B_u(i))...
    %             + T(i).*(-w).*(1/r(i)).*(r(i+1).*B_u(i+1)+2.*r(i).*B_u(i)+r(i-1).*B_u(i-1))...
    %             + T(i-1).*(w).*(1/r(i)).*(r(i).*B_u(i)+r(i-1).*B_u(i-1));
    %     end
    % end
    % 
    % %calculate new topography using implicit finite matrix
    % cu_new = bandthree(size(r,2),Ac,Bc,Cc,fc);

    %CALCULATE UPPER CHANNEL THICKNESS (more stable than implicit as of
    %8/16/23
    for i = 1:(size(cu,2))
        if (i == size(cu,2))
            cu_new(i) = cu_new(i-1);
        elseif (i == 1)
            cu_new(i) = cu(i) + (w).*(1/r(i)).*((r(i+1).*B_u(i+1)+r(i).*B_u(i)).*(T(i+1)-T(i))...
                +(r(i+1).*D_u(i+1)+r(i).*D_u(i)).*(cu(i+1)-cu(i)))...
                +(delt.*(Q*(r_o.^2./20000.^2).*exp(-r(i).^2./20000.^2)));
        else
            if (r(i) < r_o)
                cu_new(i) = cu(i) + (w).*(1/r(i)).*((r(i+1).*B_u(i+1)+r(i).*B_u(i)).*(T(i+1)-T(i))...
                    -(r(i).*B_u(i)+r(i-1).*B_u(i-1)).*(T(i)-T(i-1))...
                    +(r(i+1).*D_u(i+1)+r(i).*D_u(i)).*(cu(i+1)-cu(i))...
                    -(r(i).*D_u(i)+r(i-1).*D_u(i-1)).*(cu(i)-cu(i-1)))...
                    +(delt.*(Q*(r_o.^2./20000.^2).*exp(-r(i).^2./20000.^2)));
            else
                cu_new(i) = cu(i) + (w).*(1/r(i)).*((r(i+1).*B_u(i+1)+r(i).*B_u(i)).*(T(i+1)-T(i))...
                    -(r(i).*B_u(i)+r(i-1).*B_u(i-1)).*(T(i)-T(i-1))...
                    +(r(i+1).*D_u(i+1)+r(i).*D_u(i)).*(cu(i+1)-cu(i))...
                    -(r(i).*D_u(i)+r(i-1).*D_u(i-1)).*(cu(i)-cu(i-1)));
            end
        end

    end

    
    dT_dt = (Tnew-Told)./delt;
    dcu_dt = (cu_new-cu)./delt;
    P_l = ((rho_l-rhow).*g.*Tnew) + (rho_u-rho_l).*g.*cu;

    % Cl Calculation, Force c_i + T - cu = cl
%     c_total = cu_new+cl;
    for i = 1:(size(cl,2))
        if (cl(i) < 1) %old: 
            cl_new(i) = 10;
        else
            if (i == size(cl,2))
                cl_new(i) = cl_new(i-1);
            else
%                 cl_new(i) = a_0 + T(i) - cu(i) + delt.*E.*P_l(i);
                cl_new(i) = cl(i)+delt.*(dT_dt(i) - dcu_dt(i) + italic_f.*P_l(i));
            end
        end
    end


    
    cu = cu_new;
    cl = cl_new;


    % %Calculate volume flux between index 2 and 3 and store it for each time
    % %step
    % vel_index = 3;
    % arb1 = ((rho_l-rhow) .* ( (T(vel_index) - T(vel_index-1) )./delr) + ((rho_u-rho_l) .* (cu(vel_index) - cu(vel_index-1)) ./ delr));
    % arb2 = (rho_u-rhow) .* (T(vel_index) - T(vel_index-1) )./delr; %this is the term for derivative of pressure in upper channel dPu/dr
    % 
    % arb3 = (mu_u .* (cl(vel_index).^2)) ./ (2.*mu_u.* (cu(vel_index).*mu_l + cl(vel_index).*mu_u));
    % arb4 = (mu_l.*(cu(vel_index).^2) + 2.*mu_u.*cu(vel_index).*cl(vel_index)) ./ (2.*mu_u.* (cu(vel_index).*mu_l + cl(vel_index).*mu_u));
    % 
    % V_u(t./delt+1) = (g/mu_u).*arb2.*(cu(vel_index).^3./6)...
    %     +(g.*(-1.*arb1.*arb3 - arb2.*arb4).*(cu(vel_index).^2./2)); %velocity flux at r(2) based on topography and channel difference between 2 and 1

    %loop counter
    if (mod(t,delt*10000) == 0) %original is 10000
        t

        %20:20:800
        figure
        hold on
        for vel_index = [(150*(1/(delr/1000))) (780*(1/(delr/1000)))] %655 azores, 695 afar, 780 Iceland
            zu = 0:1:cu(vel_index);
            zl = 0:1:cl(vel_index);

            %Upper Velocity: Setup arbitrary variables (for clarity) then calculate
            arb1 = (rho_u .* ( (T(vel_index) - T(vel_index-1) )./delr) + ((rho_u-rho_l) .* (cu(vel_index) - cu(vel_index-1)) ./ delr));
            arb2 = rho_u .* (T(vel_index) - T(vel_index-1) )./delr;

            arb3 = (mu_u .* (cl(vel_index).^2)) ./ (2.*mu_u.* (cu(vel_index).*mu_l + cl(vel_index).*mu_u));
            arb4 = (mu_l.*(cu(vel_index).^2) + 2.*mu_u.*cu(vel_index).*cl(vel_index)) ./ (2.*mu_u.* (cu(vel_index).*mu_l + cl(vel_index).*mu_u));

            %velocity of upper channel in m/s
            v_u = (g/mu_u).*arb2.*(zu.^2./2)...
                +(g.*(-1.*arb1.*arb3 - arb2.*arb4).*zu);

            %Lower Velocity: Setup arbitrary variables (for clarity) then calculate
            arb3l = (mu_l .* (cu(vel_index).^2)) ./ (2.*mu_l.* (cl(vel_index).*mu_u + cu(vel_index).*mu_l));
            arb4l = (mu_u.*(cl(vel_index).^2) + 2.*mu_l.*cl(vel_index).*cu(vel_index)) ./ (2.*mu_l.* (cl(vel_index).*mu_u + cu(vel_index).*mu_l));

            %velocity of lower channel in m/s
            v_l = (g/mu_l).*arb1.*((zl - cl(vel_index)).^2 ./ 2)...
                +(g.*(arb2.*arb3l + arb1.*arb4l).*(zl-cl(vel_index)));

            vu = v_u.*speryr.*(1e3); %convert to mm/yr
            vl = v_l.*speryr.*(1e3);

            %Calculate shear stress with depth
            dvu_dz = (g/mu_u).*arb2.*(zu)...
                +(g.*(-1.*arb1.*arb3 - arb2.*arb4));
            dvl_dz = (g/mu_l).*arb1.*((zl - cl(vel_index)))...
                +(g.*(arb2.*arb3l + arb1.*arb4l));



            temprad = r(vel_index)./1000;

            plot(vu, zu./1000, 'red', 'LineWidth', 2)
            plot(vl, (zl+cu(vel_index))./1000, 'blue', 'LineWidth', 2)


            ptermu = g.*rho_u .* (((T(vel_index) - T(vel_index-1) )./delr) +(T(vel_index)./r(vel_index)));
            pterml = g.*rho_u .* (((T(vel_index) - T(vel_index-1) )./delr) +(T(vel_index)./r(vel_index)))...
                + (g.*(rho_u-rho_l).* (((cu(vel_index) - cu(vel_index-1) )./delr) +(cu(vel_index)./r(vel_index))));

            %                 dtauu_dz = gradient(tau_u);
            %                 dtaul_dz = gradient(tau_l);

            outstr = sprintf('Radius: %d, PuTerm: %d, PlTerm: %d\n',r(vel_index),ptermu,pterml);
            fprintf(outstr)

            title('Velocity-Depth Profile Plume Model')
            % xlabel('Distance + Shear Stress tau_r_z (Div 100)')
            xlabel('Distance (km) + Velocity (mm/yr)')
            ylabel('Depth (km)')
            set(gca, 'YDir','reverse')

        end
        hold off
        % tau_u = mu_u .* dvu_dz; %dz is 1 meter, no need to include
        % tau_l = mu_l .* dvl_dz;
        % %z_tau_u = diff(zu)./2;
        % %z_tau_l = diff(zl)./2;
        % % tauu = tau_u ./ 10000;%for plotting
        % % taul = tau_l ./10000;
        % %                 plot(tauu+temprad, zu, 'black', 'LineWidth', 2)
        % %                 plot(taul+temprad, zl+cu(vel_index), 'black', 'LineWidth', 2)
        % %                 plot(tau_u, zu, 'black', 'LineWidth', 2)
        % %                 plot(tau_l, zl+cu(vel_index), 'black', 'LineWidth', 2)
        %
        % dpu_dr = g.*rho_u .* (T(vel_index) - T(vel_index-1) )./delr;
        % dpl_dr = g.*(rho_u .* ( (T(vel_index) - T(vel_index-1) )./delr) + ((rho_u-rho_l) .* (cu(vel_index) - cu(vel_index-1)) ./ delr));
        % pu = rho_u .* g .* T(vel_index);
        % pl = (rho_u .* g .* T(vel_index)) + ((rho_u-rho_l).*g.*cu(vel_index));

    end



end
P_l = ((rho_l-rhow).*g.*T) + (rho_u-rho_l).*g.*cu;



